The generalized Riemann problems for hyperbolic balance laws: A 
unified formulation towards high order^ 

Jianzhen Qian^, Jiequan Li^'*, Shuanghu Wang^'^ 

"■Institute oj Applied Physics and Computational Mathematics, Beijing 100088, China 
^ School of Mathematical Science, Beijing Normal University, Beijing, 100875, P. R. China 
'^Key Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, 
fTj ■ Beijing 100088, China 

o 

(N 



(N 



Abstract 

The Generalized Riemann Problems (GRP) for nonlinear hyperbolic systems of balance laws 
in one space dimension are now well-known and can be formulated as follows: Given initial- 
data which are smooth on two sides of a discontinuity, determine the time evolution of the 
^ \ solution near the discontinuity. In particular, the GRP of {k + l)th order high-resolution is 
^ ' based on an analytical evaluation of the time derivative up to /cth order, which turns out to 
(-H ■ be dependent only on the spatial derivatives up to /cth order. While the classical Riemann 
problem serves as a primary "building block" in the construction of many numerical schemes 
(most notably the Godunov scheme), the analytic study of GRP will lead to an array of 
"GRP schemes", which extend the Godunov scheme. Currently there are extensive studies 
on the second-order GRP scheme, which proves to be robust and is capable of resolving com- 
plex multidimensional fluid dynamic problems [M. Ben-Artzi and J. Falcovitz, "Generalized 
^ . Riemann Problems in Computational Fluid Dynamics", Cambridge University Press, 2003]. 
OS ! More general formulation of the second-order GRP solver can be found in [Numer. Math. 
^ \ (2007) 106:369-425], but still confined with a class of "weakly coupled systems". In this pa- 
2 ' per, we provide a unified approach for solving the GRP in the general context of hyperbolic 
cr^ ■ balance laws, without weakly coupled constraint, towards high order accuracy. The derivation 
of the second-order GRP solver is more concise compared to those in previous works and the 
third-order GRP (or quadratic GRP) is resolved for the first time. The latter is shown to 
^ \ be necessary through numerical experiments with strong discontinuities. Our method relies 
heavily on the new treatment of the rarefaction wave. Indeed, as a main technical step, 
the "propagation of singularities" argument for the rarefaction fan, is simplified by deriving 
the L(Q)-equations, an ODE system for the "evolution" of the "characteristic derivatives" 
in x-t space for generalized Riemann invariants, with aid of the generalized characteristic 
coordinates. The case of a sonic point is incorporated into a general treatment. The accu- 
racy of the derived GRP solvers are justified and numerical examples are presented for the 
performance of the resulting scheme. 

Keywords: Generalized Riemann problem. Hyperbolic balance laws, GRP solver, 
Riemann invariants 
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1. Introduction 



In this paper we consider the generalized Riemann problem (GRP) for hyperbohc balance 
laws 

where U = (ui, ■ ■ ■ , n^) is the unknown variable with F = (/i,--- , fm) being the flux 
functions, and H{x, U) is a source term resulting from geometrical or physical effects, x is 
the spatial variable and t is the time variable. In this study, we will concentrate on the 
numerical aspect of (11.11) . rather than important theoretical issues such as well-posedness 
and solution structures. 

In the development of numerical techniques approximating solutions of fll.ll) . the finite 
volume scheme plays absolutely indispensable role, wherein one of most crucial ingredients 
is the construction of numerical fluxes and it boils down to the resolution of associated 
(generalized) Riemann problems at each computational cell interface. Specifically, we denote 
by Ij = [xj_i/2, Xj_|_i/2], Ax = Xj+i/2 — the computational cell numbered j, and by 

{tn}'^=o the sequence of discretized time levels. At = tn+i — tn- The finite volume scheme is 
then constructed by integrating the governing equations (11.11) both in space and time over 
the control volume Ij x [tn+i,tn], yielding 

= U-- - + ^tHl (1.2) 



where 



U^ = ^ Uix,tn)dx (1.3) 

is the average of U{x^t'^) over the cell Ij. The remaining terms in (II. 2p . -^j!^i/2 ^^"^ -^7 ' 
are the temporal average of F{U{x,t)) along the interface x = and the space-time 

integral average of the source II{x, U), i.e., 

F?+i/2 = ^J F{U{x^^.,t))dt, (1.4) 

"'-MA^l I , (1-5) 

A numerical scheme is obtained if one can supply suitable approximation for 
Hj with given data U{x, tn) at t = Formally, for the Godunov-type schemes, this usually 
consists the following three procedures. 
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A Data reconstruction: Based on the cell average values t/j^, reconstruct the initial data 
U{x,tn) as piece- wise smooth distribution, being constant or polynomial in each cell 

/,. 

B Solution evolution: Solve the (generalized) Riemann problem at each cell interface 
X = to evolve the solution. 

C Numerical approximation: Take the numerical integration in (11 ■4p and (II. 5p to get 
under suitable CFL condition. 

For the notable Godunov scheme Il0| and higher order schemes using Riemann solvers 
such as MUSCL 2J, |25| and TVD [ll|] schemes, the classical Riemann problem is solved 
in each cell interface to evolve the solution from to t""*"^. The corresponding initial data 
are taken as interface limit values of U{tn,x) in the neighboring cells. As an extension, the 
GRP scheme assumes the piecewise smooth initial data and evolves solutions by analytically 
solving the generalized Riemann problem at each cell interface with at least second order 
accuracy. Currently, the second-order GRP scheme has already been exploited and put into 
use for several compressible fluid models [H, [i], [^, 0, M, 0, [l6|, 27, 28|. Let us now outline a 
standard process for its implementation. Assume the data at time t = t"" is piece-wise linear 
with slope 0"", i.e. on Ij we have 

U{x, n = + aj{x - X,), X e (x,_i/2, x,+i/2). (1.6) 

Then the Godnuov-type scheme of second order takes the form 



rrn+l _ rrn ( T^n+1/2 T^n+l/2\ / n+1/2 r7-n+l/2\ 



where the following notations are used 

f;:i12 = nuTJ/i). ki'i = f^^i/?)' (i-^) 

and U^^lj^ is the mid-point value or the average of ?7(xj+i/2,i) over time interval [tn,t„+i]. 
For simplicity, the source term is currently discretized with an interface method, which is 
the trapezoidal rule in space and the mid-point rule in time 0, 13| in order to keep second 
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order accuracy. The central issue is how to obtain the mid-point value Uj^^j^, which is 
formally approximated by the Taylor expansion (ignoring the higher order terms) 



2 \dt 



-1/2 



where 



The value f^jYi/2 obtained by solving the associated Riemann problem for the homogeneous 
hyperbolic conservation laws as used in the (first order) Godunov scheme [l^. The main 
ingredient lies upon the calculation of the instantaneous time derivative (^)j_|_i/2- Even 
in the Godunov scheme, the time derivative {^)^_^_i/2 should be properly treated once the 
source term is present, which makes the solution evolve non- uniformly. 

For the solution U being smooth near the grid point (xj+i/2,tn), it follows directly from 
flLTjl that 

(f (f ^(^-'^•^-"^)' 

However, for the generalized Riemann problem including singularity at grid point {xj^i/2, tn), 
fll.lip is no longer valid, even for scalar cases, because there exists nonlinear waves (rarefac- 
tion waves or discontinuities) issuing from the singularity point (xj+1/2, tn)- Indeed, thinking 
of the initial data (11. 6p with non-zero slopes as a perturbation of piecewise constant Rie- 
mann initial data and the source term S{x, U) as a perturbation of the homogenous system 
of equations, the GRP solution is a perturbation of that of the associated Riemann problem 
at least in the neighborhood of the singularity point. It turns out that the GRP solution 
consists of, for a short time following the "disintegration" of initial discontinuity, the curvi- 
linear rarefaction wave and the discontinuities (contact discontinuity or shock wave) with 
time varying speed 0, Chap. 5]. 

The solution U together with its derivatives may undergo a jump discontinuity across 
each wave. Hence, in order to solve the generalized Riemann problem, it requires one to 
explore the mode of the discontinuity for the derivatives coming along with each wave, 
which is in fact described by a set of linear algebraic equations. This bears an analogy to 
the resolution of classical Riemann problem, which involves exploring the relation, usually 
described by a one parameter curve, between the two states of U connected by each wave. 
Indeed, the treatment for capturing the discontinuities of the derivatives across each type 
of waves can be sketched out as follows. 

A Since the generalized Riemann invariants are transported in the transversal direction 
of the rarefaction fan, it is natural to use them for studying the variation of the deriva- 
tives across the rarefaction wave. Actually, the directional (emanating characteristic 
direction) derivatives of the generalized Riemann invariants are determined by their 
values on either side of the wave. 
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B The generalized Riemann invariants, which remain continuous across corresponding 
contact discontinuities, are differentiated in the direction of the discontinuity (charac- 
teristic). 

C For the shock wave, the identities imphed by the Rankine-Hugoniot conditions are 
differentiated along the shock trajectories. 

As indicated in the previous works, the most technical step lies on the treatment for rar- 
efaction fan, which relies on the analysis in term of "characteristic coordinates" . 

The methodology for resolving the generalized Riemann problem is originated in [H, Q, Isf, 
wherein the original GRP is designed for the compressible fluid flows with two related 
Lagrangian and Eulerian versions. See also the recent textbook j3] for detailed discussions. 
The Eulerian version is always derived by using the Lagrangian case. The transformation is 
quite delicate, particularly for sonic cases, because it becomes singular at sonic points. The 
direct Eulerian version, more flexible for applications, is developed recently in the context 



of shallow water equations [16|, planar compressible fluid and the compressible fluid flows 
[l, The approach for solving GRP therein, being ready to handle any strict hyperbolic 
system endowed with a coordinate system of Riemann invariants (in particular, the two 
equations system), is extended to handle more general weakly coupled systems (in the sense 
of (gI, Def. 21]) having only a "partial set" of Riemann invariants. The common point of the 
above systems is that the generalized Riemann invariants (GRI) are coupled in a manner 
that enables a "diagonalized" treatment. Although many physical systems, including the 
compressible fluids flow system, belong to such a class of systems, the existing methods for 
deriving a second-order GRP solver are rather complicated, which prevents it from practical 
use in many ways. For example, the treatment of rarefaction relies heavily on the explicit 
formulation of the Asymptotic Characteristic Coordinate (ACC), which depends on the EOS 
(equation of state) of the fluid in turn and is sometimes hard to derive. Besides the ACC 
is not easy to be written out explicitly for higher order GRP solvers that are particularly 



useful in capturing the propagation of entropy wave 20| (see also Fig. 18. 4p . Other closely 
related efforts can be found in 0, [l^ using the approach of asymptotic analysis for the 
resolution of generalized Riemann problems, and in jlil . [iil . [sj (and the references therein) 
for approximate Godunov-type high order solvers. The solvers in [2^, @] corresponds to 
the acoustic case and they fails for resolving strong discontinuities, even with very high 
order accuracy (this point is confirmed through a numerical experiment Fig. I7.4p . Hence 
it is absolutely necessary to develop the high order (at least third order) GRP scheme by 
resolving nonlinear wave patterns each computational grid point analytically, in addition to 
provide an acoustic approximation as the jump there is weak. 

Therefore we present a unified approach in this paper, still direct Eulerian, to resolve the 
GRP for general systems of hyperbolic balance laws fll.ip . The weakly coupled constraint 
is not required here. The solver for second-order GRP (linear GRP) as well as third-order 
(quadratic GRP) are derived. This paper provides a simplified treatment for the main step, 
resolution of the rarefaction fan. Indeed, it is carried out by first deriving the system of 
transport equations for the generalized Riemann invariants and then deriving the "evolution" 
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equations, labeled as the L(Q)-equations, for their characteristic derivatives in x-t space 
with aid of the generalized characteristic coordinate (GCC). This is based on the following 
observations. These characteristic derivatives, and hence the resulting "evolution" equations 
for them, are independent of the choice of the auxiliary GCC. Thus the explicit expression of 
the GCC is not required. More importantly, only in the emanating characteristic direction do 
the derivatives of the GRI (the solution U) of any order exist and remain continuous across 
rarefaction fans. No additional assumption is required for the regularity of the solution U. 
Indeed, the above observations are the reasons why we can derive the third-order (or higher 
order) GRP solver without many difficulties. By referring to Section H] for the resolution 
of the contact discontinuity and the shock wave, the spatial derivatives of the solution, 
from which the instantaneous time derivatives follow, are obtained by solving a simple 
system of linear equations in the intermediate regions of the waves. The case of sonic 
point is handled by supplementing an additional freedom using the differential relation of 
U along the emanating characteristic direction. A special case frequently occur during the 
numerical application of the GRP scheme is the acoustic case: the initial values of U are 
continuous at the singularity point. This case is comparatively easy to handle and requires 
less computation cost. 

Although this paper focus on exploring solvers for the second-order linear GRP and the 
third-order quadratic GRP, higher order GRP solvers can be derived with the same method- 
ology and a multidimensional extension can be pursued in a forthcoming work The 
resulting GRP solvers consist of two steps: (i) the classical Riemann solver; (ii) calcula- 
tion of instantaneous time derivatives of U. As indicated by the solvers. Step (ii) can be 
straightforward once the full Riemann solution is obtained. Besides, in both steps, only the 
limiting values of U and its spatial derivatives at two side of the singularity are used, and the 
resulting linear (resp. quadratic) GRP solver leads to second (resp. third) order accuracy 
in time approximation to U regardless its initial distribution. 

This paper is arranged as follows. In Section 2, a basic setup for the system and the 
GRP are presented. The resolution of rarefaction wave and discontinuity waves, including 
the contact discontinuity and shock wave, are detailed in Sections 3 and 4, respectively. We 
conclude the resolution of GRP in Section 5 and the acoustic approximation in Section 6. 
As an application example, in Section 7, we derive the GRP solvers for compressible variable 
duct flow system and show the solvers' accuracy by several tests. Finally, in Section 8, the 
GRP solvers are used to construct one-step high order numerical scheme and a few 1-D 
numerical test cases are presented. 



2. Basic setup for the system and the GRP 

As a basic setup, we assume fll.ip is hyperbolic in the sense that the Jacobian A[U) = 
^^Jy^ of F{U) has m eigenvalues 



Ai < As < ■ ■ ■ < A 




6 



The set of left (right) eigenvectors Lk (Rk) (associated with X^, k = 1, ■ ■ ■ ,m) are hnearly 
independent. The kth characteristic field Xk can be either genuinely nonlinear in the sense 
of Vi/Afc ■ i?fc 7^ 0, or linearly degenerate Vc/Afc ■ Rk = 0. 

Now let us state the generalized Riemann problem. It is defined as the initial-value 
problem for system fll.ip . subject to the initial data 

where P±{x) are vectors, whose components are the smooth functions. As illustrated in 
Section [H the initial structure of the solution is determined by the associated Riemann 
problem: 

at + & u, ^^^^^ 

f/^(x,0) = f/±, ±x>0, 

where U± are the limiting values of P±{x) at x = 0, i.e. U± = -P±(0^). We call the solution 
of (12. 3p the associated Riemann solution of (II. ip and (12. 2p . 

Assumption 2.1. The Riemann problem 112. 3\) is uniquely solvable, and the solution to 
\2.3\) consists of m waves Fi, ■ ■ ■ , r^. The wave Tkfl < k < m) is an admissible shock, 
a contact discontinuity, or a rarefaction wave associated with the kth characteristic field Xk- 

Note that the above assumption does not mean we are confined with strict hyperbolic 
systems that endowed with distinct eigenvalues. 

Denote by R'^{x/t, f/_, f/+) the Riemann solution of (12.31) . Then we have the following 
proposition. 

Proposition 2.1. Let U{x,t) be the solution to the generalized Riemann problem U.l\} and 
112. 2\} . Then for every fixed direction 6 = x/t, 

\imU(et,t) = R^(e,U^,U+). (2.4) 

This implies that the wave configuration for the generalized Riemann problem U.l\) and Ii2.2\} 
is the same as that for the associated Riemann problem 112.!^) around the singularity {x, t) = 
(0,0+). 

Proposition l2.1l is illustrated schematically in Fig J2.1[ The solution of (12. 3 p is self-similar, 
and hence the waves are centered. Correspondingly, the waves for (II. ip are curved (See [3] 
for more detailed descriptions). 

We emphasize that the solution U is smooth in the intermediate regions of these waves 
and along each emanating characteristic curve in the rarefaction fan (up to the singularity 
(0,0"*")). To approximate U along t-axis with kth order accuracy, we can use the Taylor 
expansion 

U{x = 0, t) = f/(0, 0+) + E ^^(0' + ^(^'^')- (2-5) 




U-(x,t) 



U+(x,t) 



(a) 



(b) 



Figure 2.1: Wave configurations: (a) Wave patterns for tlie GRP with initial data U{0,x) = P-{x) for 
X < and U{0, x) = P+{x) for x > 0, U± = P±(0±). (b) Wave patte rns for the associated Riemann problem. 



A solver of the GRP is actually that of evaluating the instantaneous time derivatives 



-(0,0" 



lim-r^(0,t), 



t > 0. 



(2.6) 



For convenience, we label the problem of evaluating 02. 6p with i = 1 (resp. £ = 1,2) as the 
linear GRP (resp. quadratic GRP), or LGRP (reps. QGRP) for short. As mentioned in the 
introduction, this paper concentrates on QGRP. 



3. Resolution of curved rarefaction waves 

As pointed out earlier, the main feature of the GRP is the resolution of rarefaction 
waves and the main ingredients are the Riemann invariants and characteristic coordinates. 
Let us consider to first derive the set of transport equations for the (generalized) Riemann 
invariants in a general setting. For this purpose, we rewrite (11. 1[) as a nonconservative form 



by recalling A{U) = dF{U)/dU. Muhiplying (E]) by L = (Li, 
follows that 

dU dU 

■ ■ , Am)- If there exists a set of variables w 
m), then (13. 2 p is equivalent to 

dw . dw 



where A = diag(Ai, 



(3.1) 

Lrn) from the left, it 

(3.2) 
satisfying 



,Wr. 



dt ' dx 



LH{x,U). 



(3.3) 



Indeed, w is a complete set of Riemann invariants. Unfortunately, most of the systems (II. ip 
with m > 2, including the full system of compressible Euler equations, do not admit such 
a set of Riemann invariants. We thus turn to exploring the generalized Riemann invariants 
(GRI). 
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3.1. The generalized Riemann invariants (GRI) 

Let w = {wi, ■ ■ ■ , Wm-i) be the generalized Riemann invariants of the kth characteristic 
field. By recalling the definition of GRI 22| , we have 

VuWi-Rk = 0, i = l,---,m-l. 

Hence, there exists an invertible (m — 1) x (m — 1) matrix K, such that 

KVuW = (Li, ■ ■ ■ , Lfc_i, Lfc+i, ■ ■ ■ , L^)^ =: L'-'\ (3.4) 

Multiply (13. ip by K'Vuw from the left yields the following proposition. 

Proposition 3.1. Let w = [wi, ■ ■ ■ , Wm-i) be the GRI of the kth characteristic field. Then 
in any smooth region of U there holds 

^ + B^''\U)^ = L^''^H{x,U), (3.5) 

where 

B^^\U) = K-'A^'^^K, A(^) = dzagiX,, • • ■ , A^.i, A,+i, ■ ■ ■ , A^)^, (3.6) 
and K is determined by [3^. 



Roughly speaking, Proposition 13.11 implies that the generalized Riemann invariants of the 
kth characteristic field are transported along the direction different from A^. The following 
useful corollary is straightforward from Proposition 13 . 1 1 for the resolution of rarefaction wave. 
See Remark 13.11 (ii) below. 

Corollary 3.1. Let Tk be a characteristic curve associated with A^. If U is continuous 
and piecewise smooth with Tk being a weak discontinuity curve, then dw/dt and dw/dx in 
Proposition \3.1\ remain continuous across Tk. 

3.2. Generalized characteristic coordinates (GGG) 

As mentioned in the introduction, the characteristic coordinates, defined as the integral 
curves of the characteristic equations, play an important role in the resolution of rarefaction 
waves. In the region of a rarefaction fan, they work similarly to the usual polar coordinates 
to single out singularities. 

Assume that is a rarefaction wave associated with A^ and denote by ^/^(x, t) (resp. 
Uji{x,t)) the state U on its left (resp. right) side. See Fig. 13.11 To simplify notations, we 
write A and B below for A^ and B^''^ in Proposition 13. Ij respectively, for a fixed k. 

Let C~: (3{x,t) = (3 and C+: a{x,t) = a, (3 E [(3l,(3r], — oo < a < be the integral 
curves of the following equations, respectively, 

dx dx 

Here, different from the previous works [l|, 0, 3, Q, S 0], /i in (13. 7p is not required to be 
an eigenvalue of A{U). In fact, it can be the slope of any family of transversal curves 
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C- :!3 = !3l 




UR{x,t) 



Figure 3.1: Generalized characteristic coordinates in a rarefaction fan Ffc. 



different from A. For example, /i = —K Tlie GCC used here is a relaxed version of the 
afore-mentioned ACC. Moreover, /3 and a are denoted as follows: /3 is the initial value of 
the slope A at the singularity (x, t) = (0, 0) and a for the transversal characteristic curves 
is the x-coordinates of the intersection point with the leading /3-curve: /3 = Pl- 

The coordinates {x, t) in the "triangle" sector of the centered rarefaction wave shown in 
Fig. 13.11 can be expressed in terms of a and (3, 



which satisfy 



Denote 



Then we have 



X = x{a, (3), 



dx ^ dt 
da da ' 



— A— 

dt dx' 



t = t{a, /3), 

dx dt 
dp^^dp' 

^ d d 



d_ 

da 



dt 



In particular, as a = 0, we have 



da 
dt 



Dy 



d 



dt 
dt 



dx 



d(3 d(3^''' 



(3.8) 
(3.9) 
(3.10) 
(3.11) 

Here we remind that, as a basic assumption in Section |2l the solution U is smooth along 
each characteristic curve C~ : ^ = P inside a rarefaction fan up to the singularity and D^xU, 
for any i > 1, takes finite value at a = 0. 

It follows, by differentiating the first equation in f l3.9p with respect to f3, the second with 
respect to a and then subtracting, that the function t = t{a, /3) satisfies 

dH dX dt dfi dt 



(0,/3) = 0, f3L<(3<(3_ 



(/x-A) 



dad/3 df3 da da df3 ' 



(3.13) 
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Setting a = and using fl3.12p . one obtains 



d_ 



dt 



(0,(3). 



We continue to make differentiation of f l3.13p witli respect to a to obtain 

dH d , , , dH d^X dt dX dH 



(3.14) 



(/i-A) 



da^dp 



A(^_A)^!^ + + 

da dad 13 dad (3 da d(3 da^ 

d'^li dt dfi d'^t 



da^ d(3 da dad (3 
Recalling f l3.12p and fl3.14p as well as noticing 



(3J5) 



9^ A = 4r5«.A 



dad/3 d/3 \da 
we take a = to obtain 

dH 



dH dt d 



d_ 

dp 



da' 



r(0./?) 



/i — A da'^ 



ji — X \da 



dt 



(0,/3) 



+ 



d 



fx — X d/3 



(^aA) 



dt_ 

da 



(0,/3). 



(3.16) 



The equations 03.141) and (13.160 . for dt/da{0, /3) and d'^t/da'^{0, /3) repectively, are crucial 
for deriving the L((5)-equations in next subsection. 

3.3. The L{Q) - equations 

In this subsection, we shall derive the linear differential ordinary equations for DxwlO, /3) 
and D\w{Q,/3) with respect to /3, namely, the L-equations and Q-equations, respectively, of 
the GRI. Precisely, we have the following proposition. 

Proposition 3.2 (L((5)-equations). The w in Proposition \3. 1\ satisfies the L-equations: 

d 



d[3 



[Z^A^O, (3)] = {XI - B)-' {Dxw - L^^'^H) , 



(3.17) 



and the Q-equations: 



d 



— [DXO, (3)] =2(A/ - B)-^Dlw + 2Dx{XI - By^D^w - 2D, [{XI - By^L^'^^H] 



+ ^{DxX) [{XI - B)-'D,w - {XI - B)-'L^'^H] , 



(3.18) 



for/3e[/3L,/3j,]. 
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Remark 3.1. (i) Note that jS.ll^ and j<§]) are the linear ordinary differential equations 
for Dxw and D\w, respectively. Equivalently, by integration, we can formulate ( [5'. i7p as 

DM0, f3) = C^''^DM0,f3L)+SC^'\ (3.19) 

and l{3.18\} as 

DlwiO, (3) = Q^'^Dlw{0, (3l) + SQ^''\ (3.20) 

where C^^^ and Q^^^ are both the (m — 1) x (m — 1) matrices andSC^'^\ SQ^^^ are the (m — 1) 
vectors. 

In particular, for most physical systems, the S in ( (J.y/[ ) is a sub-triangular matrix, so 
are B and (A - B)"^ in and ^EHE)- Hence \3.1!J\) and l[3.2U\) can be obtained by 

integrating component by component of D\w in \3.11^ and D\'w in ^3.18\] . respectively. 



(a) Corollary \3.1\ ensures that d^w remains continuous across both the head P-curve: 
(3 = Pl and the tail (3-curve: (3 = /3r. Moreover, \3.1'1\) is equivalent to the following 
equation of d^w 



^ [dMO, /?)] = (A/ - Br [-^Bd.^ - ^(L^'^H 



which can be formulated as 



(3.21) 



(3.22) 



However, we can not derive an equation analogous to \3.21]) for dl'w(0,l3) since (9^iw(0,/3) 
for (3 G (/3l, I3r) does not take a finite value in general. 



Proof of Proposition 13. 2L We make use of the regularity of Riemann invariant w. Let 
us first differentiate w witli respect to a and /3 to get 



d dw d f dt ^ 



da 8(3 da \d(3 " 

Similarly, one has 

d dw d f dt ^ 

— I —Dxw 



3(3 da 3(3 \da 
Subtracting these two equations yields 



dt d ^ dt d ^ 
Using fl3.12p and fl3.14p . one can obtain 
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dH ^ dt d ^ 



dj3da ^ da 8/3 ^ 



dH 

dad [3 



{Df, - Dx) w = d^w. 



(3.23) 



(3.24) 



(3.25) 



Recall Proposition 13.11 Then we arrive at fl3.17p . 

We proceed by differentiating f l3.23p and f l3.24p with respect to a to obtain 



fdt 



^ ^ dH d ^ dt 92 ^ 



da^df} ^ dad (3 da 



3(5 da 



2 ' 



and 



^3 



f dt 



d(3da' 



-w 



d(3da \da 
dH 



d f dH 



dt d 



dp Kda'^^^'^ ~^ da da^^"^ 



d/3da' 



dH d 
d^d^ 



DyW + 2 



dH dt 

dad 13 da 



r.2 fdtV d ^2 

Diw + — —Diw 



\da J dp 



Subtract the above two equations and then set a = to yield (using ( I3.12p and fl3.14p again) 

d'^t d'^t d 



dp da 



da'^ dp 



Recalling fl3.25p . it follows that 
f dtV d 



- 



\da J dp 



2 

+ 



- D^{d.,w) 



dt\ 

da J 

dH , dH 
(/i- A) 



f dt 



- 



dpda^ 



/i — A \da 
dxW. 



da 



Inserting f l3.16p into the last term of the above equation, we can obtain (after suitable 
reduction) 



— [DX0,/3)] = 2Dy^{d^w) + —{Dx\)dxW. 
Recalling Proposition 13. II once more, we obtain f l3.18p . 



(3.26) 
□ 



4. Resolution of curved discontinuities 

In this section, we resolve the curved discontinuity wave, which can be a contact dis- 
continuity or a shock wave. Let Vk be the discontinuity wave and denote by UL{x,t) (resp. 
UR{x,t)) the state U on its left (resp. right) side. 
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4.1. The contact discontinuity 

Assume for the present that the Tk is a curved contact discontinuity. We use the same 
notations as in the previous section. The propagation speed of is A by suppressing the 
subscript for simphcity in notations. 

A significant feature of contact discontinuity is that the generahzed Riemann invariant 
remains continuous across the wave. Thus we take differentiations of w along the trajectory 
of Tk to obtain 

Di{w{Un))=Di{wm), (4.1) 

for ^= 1,2. 

By recalhng (13. 5p . we have 

Dx{wiU)) = [(A/ - B)Wuw]idM) + L^'^H. (4.2) 

Thus, while i = 1, (14. ip is equivalent to 

[{XI - B)Vuw]R{d,U)n - [(AJ - B)Vuw]L{d,U)L = -{L^''^H)r + {L'^'^H)l. (4.3) 

4-2. The shock wave 

Now, let us assume Tk is a curved shock wave with propagation speed denoted by cr. 
Then along the shock trajectory, the Rankine-Hugoniot relation reads 

F{Ur) - F{Ul) = a{UR - Ul). (4.4) 

Denote Do- = ^+^^- By taking the directional derivative of (14. 4p along the shock trajectory 
Ffc, one can get 

Di{F{UR) - aUn) = Di{F{UL) - aU^), (4.5) 

for 1,2. 

While £ = 1, by noting that 

D„{F{U) - aU) = {A- aI)DM - D^aU 

= -{A- alf{d^U) -{A- aI){L^''^H) - D^aU, (4.6) 

(14. 5 p is equivalent to 

- (Ar - aI)\d,U)R + {Al - aI)\d,U)L - D,a{Un - Ul) 

= -{An - aI)iL<^'^H)n + {A^ - aI){L^''>H)L. (4.7) 

An alternative approach for resolving the shock wave is by using the m — 1 Rankine- 
Hugoniot relations in the form 

*(f/L, t/i?) = 0, * = ■ ■ ■ , D), (4.8) 

which is equivalent to (14.41) . By differentiating fl4.8p along the shock wave, the relation 
equations for (D^[/)/j and (-D^[/)i(/c = 1,2) can be obtained directly. This later approach 
is usually more efficient for practical use, since a does not appear in fl4.8p now. However, 
for a general purpose, we shall use the former approach in the following discussion. 
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5. The GRP solvers 



In this section, we wiU present the fuU solver for the linear GRP and quadratic GRP. 
Indeed, since the solution U is smooth in the region on the left (resp. right) of Pi (resp. 
Pm), the time derivatives of U are thus determined by (11. ip and the initial data (12. 2p . 

For the nonsonic case, it suffices for us to determine the spatial derivatives d^U and d'^U 
in the intermediate regions of P^ {k = 1, ■ ■ ■ ,m), since the times derivatives dtU and d^U 
follows directly from (II. ip . For the sonic case that t-axis lies inside the rarefaction fan, we 
need to give an independent treatment. 

5.1. The nonsonic case 

The nonsonic case refers to the case that the t-axis is located in the intermediate regions 
of Pfc (/c = 1, ■ ■ ■ ,m). The m waves P^, k = 1, ■ ■ ■ ,m separates the half space t > into m+1 
regions. The region on the left (right) of P^ is labeled as flk-1/2 (^fc+i/2)- The associated 
state of U in Qk-1/2 is labeled as Uk-i/2- The same notation apply for the derivatives of ?7, 
such as {dxU)k-i/2- 

Now let us summarize the resolution of the linear GRP for the nonsonic case in the 
following proposition. 

Proposition 5.1 (Linear GRP: Nonsonic case). Assume that the solution of problem 
M.l\) and 112. consists of m waves Tk, k = 1, ■ ■ ■ ,m. Then the (m — 1) x m unknowns 
{dxU)k-i/2, (k = 2, ■ ■ ■ ,m) in the intermediate regions ofT^ and the number D^^ak are 
determined by the following linear algebraic system 

' {Vuw)k+i/2{d,U)k+i/2 - M^'K^uw)k-i/2{d,U)k-i/2 = SM'^''\ 
ifVk is a rarefaction wave; 

[(AfcJ - B^^^)Vuw\k+i/2{dxU)k+i/2 - [(AfcJ - B^^^)Vuw\k^i/2{d.,U)k-i/2 
= -(LWi7)fe+i/2 + (L('=)i/)fc_i/2, 

< 

ifTk is a contact discontinuity wave; 

-{Ak+1/2 - (Xfe/)2(9^t/)fc+i/2 + (Afc„i/2 - Crfc/)^(9^f/)fc_l/2 - Dak(^kiUk+l/2 " f/fc-l/2) 

= -{Ak+i/2 - afc/)(LWi7)fc+i/2 + {Ak.i/2 - ak){L(''^H)k-i/2 
if Pfc is a shock wave. 

(5.1) 

Here, the and SM^''^ are as m (EJD- (<9^f/)i/2 = {dxU^L, (9xt/)m+i/2 = {dxU)R and 

the Uk+ii2 in the coefficients are determined by i?"^(6', f/_, [/+). Having solved {dxU)k-i/2, 
the time derivatives {dtU)k-i/2 (k = 1, - ■ ■ ,m + 1) are determined by 

{dtU)k^i/2 = ~A{Uk^,/2){dxU)k^i,2 + H{x, f/fc„i/2). (5.2) 
Remark 5.1. To solve ( 15. ij) . we suggest that the unknowns be ordered as 

(■ ■ ■ ' ^a-k'^k, Uk-1/2, Uk+1/2, ■ ■ ■ ) 

if Tk is a shock wave, and use Gauss-Jordan elimination with rows partial pivoting. 
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Proof of Proposition 15.11 As illustrated previously, the solution U of (11. ip and (12. 2p 

is smooth in the regions Qk-i/2, k = 1, ■ ■ ■ ,771 + 1. In the regions Q1/2 and Qrn+1/2, the 
spatial derivatives dxU are determined by the initial data {dxU)L and {dxU)^, respectively. 
As indicated by the resolution of rarefaction wave and discontinuous waves in Sections [3] and 
m the relations between {dxU)k+i/2 and {dxU)k^i/2 are described by a set of linear algebraic 
equations. (15. ip is obtained by combining (I3.22p . (14. 3 p and (14. 7p and (15. 2p follows directly 
from dal]). □ 
To present the quadratic GRP solver, we need to give a few formulations. In the regions 
where the flow is smooth, by applying and dt to (II. ip . we have 

dtidxU) = -AdlU - dxAdxU + dxH, (5.3) 
d'^^U = -AdtidxU) - dtAdxU + dtH. (5.4) 

Inserting (II. ip and (15. 3 p into (15. 4p . dfU can be expressed as a function of t/, dJJ, dIU: 

d^U = AQ{U,dxU,dlU). (5.5) 

For the GRI w of kth characteristic fields, by noticing that 

9> = dx{Vuw)dxU + VuwdlU 

Dx.idxw) = (AfeJ - B^'^)dlw - dxB^'^dxW + dxiL^'^H), 

Dlw = D^^XJ - B^^^)dxW + (AfcJ - B^''^)Dx^{dxW) + Dx^{L<^''^H), 

we can get 

Dlw = Ml'\U)dlU + i?f (t/, dxU), 

with 

Ml^\U) = {XuI-B^^^fVuw, 

B^J^\U,dxU) = [D.^iXJ - B^^^) - (A, - B^''^)dxB^''^]dxW (5-6) 
+ (A,/ - B^'^Ydx{Vuw)dxU + (AfcJ - B^'^)dx{L^''^H) + D^^{L^^'^H). 
To resolve the shock wave, we shall use 

Dl{F{U) - aU) = D,{{A - aI)D„U - D^aU) 

= D^{-{A- alfdxU +{A- aI)H - D^aU) 
= M,{U, a)dlU - DlaU + 5,(f/, d^U, a, D„a), 

with 

M,{U,a) = {A-alf, 

Bs{U, dxU, a, D^a) = {D^A - 2D^aI)D^U + {A - aI)\dxAdxU - d^H) (5-7) 
+ {A- aI)[-DM - <yI)dxU + D^H]. 

Similar to Proposition 15.11 by combining (13.181) and (14.11) . (14.51) with ^ = 2, we have the 
following proposition for the quadratic GRP solver in nonsonic case. 
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Proposition 5.2 (Quadratic GRP: Nonsonic case). Assume that the solution of prob- 
lem U.l\} and 112. ^)) consists of m waves Tk, k = 1, - ■ ■ , m. Then the (m — 1) x m unknowns 
{d1U)k~i/2 (k = 2,--- ,m) in the intermediate regions of Tk and the number D'^^.ak are 
determined by the following linear algebraic system 



k-l/2 

T^{k),TT ICS TT^ 

k-l/2) 

if Ffc is a rarefaction wave; 



Mr{diu)k+i/2 - Mr{Uk-i/2mu)k-i/2 = 

-#([/fc+l/2, (9xt/)fc+l/2) + #(f/fc-l/2, (5xt/)fc-l/2), 

if Ffc is a contact discontinuity wave; 

M,{Uu+i/2.a){dlU)u+i/2 - Ms{Uk-i/2,(T){dlU)k-i/2 - Dl^au{Uk+i/2 - Uu-1/2) = 
—BsiUk+i/2, idxU)k+i/2, 0-k, Dfj^ak) + Bs{Uk+i/2, {.dxU)k+i/2, o"fc, D^j^ak), 

if Vk is a shock wave. 

(5.8) 

Here, the Mr^.''\u), Bi^\u,dxU), Ms{U,a) and BsiU, d^U, a, D^a) are as m (EJj and 
{dlU)i/2 = {dlU)L and {d'^U)m.+i/2 = {dlU)R. The Uk+1/2 and {dxU)k+i/2 in the coefficients 
are determined by Proposition \5.1[ Having solved {d1U)k-i/2, the time derivatives {dfU)k-i/2 
(k = I,-- - ,m + 1) can be obtained by using 

{d^U)k-i/2 = AQ{Uk-i/2, {dxU)k-i/2, {dlU)k-i/2). 

5.2. The sonic case 

As far as the sonic case is concerned, the t-axis is located inside the rarefaction wave, Tk 
for instance, and is in fact tangential to the Afc-characteristic curve. Thus, for this case, we 
need to solve Dx^U and D\JJ. Moreover, the exphcit expression of D\^U (with respect to 

is required when solving the Q-equations (13.181) . as is the main step of the QGRP solver. 

Although Dx^w and D\^w are readily obtained from (13.191) and (13.201) . we still need to 
make up an additional freedom. In fact, (13.21) implies the differential relations of U along 
the Xk characteristic curve 

LkDx,U = LkH. (5.9) 

Combining (15.91) with 

WuwDx.U = Dx,w, (5.10) 

DxJJ can be determined. 

Furthermore, applying Dx^ to (15. 9p and (I5.10p yields 

LkDlU = -Dx,LkDx,U + Dx.iLkH), (5.11) 

and 

VuwDlU = -Dx,iVuw)Dx,U + Dlw. (5.12) 
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C : f) = p. 








X 



Figure 5.1: The characteristic curve {3 ~ red hne. 



Then D\JJ can be solved by combining (15. lip and (I5.12p . 

In the sonic case, where = 0, we use the following observation for the LGRP, 



for t > 0. However, the above observation can not be used to calculate c?^?7(0, 0"*"), since 
generally neither d^U (0, 0+) nor dJJ (0, 0"'') takes finite value inside the rarefaction wave fan, 
expect for U being the GRI w. See Remark 13.11 

For the QGRP, we shall use the following method to give a second order in time approx- 
imation of U in t-axis. For any point P* = (0, At) with At being small, to evaluate U{P^), 
we need to find the initial slope /3o of the characteristic curve C which emanates from the 
singularity and goes though P^. See Fig. 15. 1[ Since for any {x(t),t) G C, 



f/t(0,0+) = D,,.t/(0,0+). 



(5.13) 



Also, by taking = 0, we have 



(5.14) 



X 



^t)= [ \ds= [ (A,,(0) + DA,Afc(0)s + O(s2))ds, 



Jo Jo 



(5.15) 



the initial slope /3* = Afc(O) can be approximated by solving 



/3, + Z}A,Afc(/3,)— = 0, 



(5.16) 



for which, we can use the Newton iteration with initial guess /3* = 0. 
Having determined U{P^) can be evaluated as 



U{P,) ^ f/(/3.) + DA,f/(/3*)At + DllJ{P,) — 



(5.17) 
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6. The acoustic approximation 

As f/_ = U+ and 9fP„(0~) 7^ 9fP+(0"^), we refer it to as the acoustic case and all waves 
Ffc are acoustic. Fixed A^, the wave degenerates to a characteristic curve and the states 
Ul, Ur on both sides of are the same. In particular, as the initial data has a small jump 
||f/„ — f/+|| ^ 1, we adopt the acoustic approximation in the sense that f/_ and f/+ are 
regarded as the same approximately. 

Now let us look at the acoustic wave r^. We use the continuity property of U and make 
differentiation along to obtain Dx^Ul = Dx^Ur. Then we proceed to use (13. ip to get 

{\I - A)L{d,U)L + Hl = {\kl - ^)i?(5.f/)R + Hr. (6.1) 

Note that Ul = Ur and recall the notation L'^'^^ = (Li, ■ ■ ■ , Lfc_i, Lfc+i, ■ ■ ■ , iii (13. 4p . 
Then we find that the (16. ip is equivalent to 

L'^^\d^U)L = L'^>'\d^U)R. (6.2) 

Moreover, applying Dx^. to (16. 2p yields 

Dx,L^'^ {{d,U)L - {d,U)R) - L(^) [{d^Ad,U)L - {d,Ad,U)R) . (6.3) 

In addition, if {dxU)^ = {dxU)R, then (16. 3p is reduced to 

L('\dlU), = L('\dlU)R. (6.4) 

Interestingly, in the course of acoustic approximation, we can obtain d^U equivalently 
by solving linear classical Riemann-type problems 

i{diU)+A{U,)i,{diU)=0, 

f,iTT( m-J ^'^-(0) ifa;<0, (6-5) 
C.^l^^Uj - I ifx>0. 

with = {U_ + f/+)/2. Note that the components of L^'^'^d^U in are nothing but the 
m — 1 generalized Riemann invariants of system (16.50 associated with A^. We also note that, 
as indicated by (16. 2p or (16. 5p . the spatial derivatives dxU are independent of the source term 
H. In addition, if dxP-{0) ~ 9^P+(0), from (16.41) . we see that d^U can also be approximated 
by solving (16. 5p with i = 2. In general, we have the following proposition. 

Proposition 6.1. For any k > 1, assume that we have 9fP_(0) = 9f.P+(0) (0 < £ < k — 1 
with d^U stands for U). Then d^U are determined by the linear system /i6.5\) with i = k. 

Remark 6.1. (i) Since U is analytical in the regions (j = 1? ■ ' ' j''^ + ly'? '^he corre- 

sponding time derivatives dfU, < i < k follow from the Cauchy-Kowalewski procedure as 
illustrated in fil /. 

(a) We note here that, under the acoustic assumption in Proposition ^ . 1\ all the approxi- 



mate DRP solvers proposed in lln . \23l . \aj are valid and are actually equivalent to the present 
acoustic GRP solvers. 
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As for the resolution of GRP fll.ip - fl2.2p . if the initial data (12. 2p has a jump discontinuity, 
we can derive the solvers analytically as in Section |5] to calculate the time derivatives of U, 
with possible acoustic approximation for a partial set of waves. This leads to the solver 
which we label as the LGRPoo (QGRPoo) solver. While the jump U+ — U^oiU is very 
small, we can use (16.51) (or possibly (16. 3p ) to calculate the space derivatives approximately. 
The resulting LGRP (QGRP) solver is labeled as the LGRPi (QGRPi) solver. 



7. An example: The variable area duct flow system 

In this section, we will take the system of variable area duct flow as an example to test 
the GRP solvers in the previous section. The flow system is 





U=\pu\, F{U)= pu'+p , GiU) = 
\eJ \{E + p)u) 

Here, p, m, e are the density, velocity and internal energy, respectively, p = p{p, e) is the 
pressure, E = p(e+ l/2u^) is the total energy and the function A{x) is the area of the duct. 
When A{x) = 1, the system (17.61) represents the planar compressible Euler equations. We 
discuss the case of polytropic gases, for which p = (7 — l)pe, where 7 is the ratio of specific 
heats. 

7.1. Formulation of the GRP Solvers 

In terms of the primitive variables Q = (p, u,p), system (17.61) can be written, for smooth 
flow, as 

Vo pc2 uj y-^pc^u) 

Here, c is the local speed of sound, given by = y. 

The system (17. 6p . or equivalently (17. 7p . possesses three eigenvalues 

A_ = n — c, Aq = li, A+ = n + c. 

The three pairs 

w^ = {S,ip), Wo = {u,p), w+ = {S,(l)) 

are the generalized Riemann invariants associated with A_, Aq, A+. Here, S = pp~'^ is the 
the entropy, and the two varibles ipy are 

2 2 

ijj = u-\ c 4> = u c. 

7 — 1 7 — 1 



20 




Figure 7.1: Typical wave configuration of tlie variable area duct flow system. 

We now start to resolve the generalized Riemann problem for (17.61) subject to initial 
data (12. 2p . Assume that the configuration is as shown in Fig |7.1l a rarefaction wave r_ 
associated with A_ moves to the left, a shock wave r+ associated with A+ moves to the 
right. For the variable U, let us denote by Ul and f//j its values on the left-hand side and 
right-hand side of the three waves, respectively. Similarly, the values of U on the two side of 
To are denoted by Ul and U^, as is illustrated by Fig. 17.11 Similar notations will be applied 
to other variables. For example, {dxU)*j^ is the value of dxU on the left-hand side domain of 
the contact discontinuity. 

As stated in Sectional to resolve rarefaction wave associated with A_, we need to use 
the associated GRI w^. Indeed, in view of Proposition 13.11 we have the following equation 
for 

By recalling 5(0,/?) = 5^, ^/'(0,/3) = 0l, A_(0,/3) = /3 in the rarefaction wave fan and 
using the L(Q)-equations in Proposition 13. 2[ we can obtain the following proposition. The 



coefficients Ai, ■ ■ ■ , Aig, Bi, ■ ■ ■ , B12, ^i(/3) and ^2(/3) are given in [Appendix A 



Proposition 7.1 (L(Q)-equations: 7 7^ 3, 5/3). Let T_ be the rarefaction wave as in Fig. 



7.1\ Then for 7 7^ 5/3, 3, we have 



7+1 



D,_SiO,l3) = A,{tljL-P)- 
Da_^(0, P) = A,{^L A.ii^L - P)^ + B,{^L - /3) + B^ii^L - Pf. 

and 

DaX0,/3) = [{^Pl - /3l)^^/^aXO,/3l) + Z^iP) - Z2{f3L)]{^L - /3)^, 
for/3e[f3L,/3l]- 
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(7.9) 



(7.10) 



Remark 7.1. (i) For the cases of'y = 5/3, 3, the L-equations of are given in Appendix A 



The corresponding Q-equations are omitted here. Besides, the formulae for dxS{0, (3) and 
9a;\l/(0, which can be obtained by integrating 113. 21\} directly, are also omitted here. 

(a) For '-f 7^ 5/3,3, the differential relation in the rarefaction wave fan associated with 

A_ 

^A_0 = -^^>A_5 + ^CU (7.11) 

7(7 - 1) 5 A{x) 

leads to 

Dx_(j) = A,{^L - + Bsii^L -13) + B.iijL - (7.12) 



By noticing A_ = "^-^il) + ^-^4>, 
we have 

^(Da_A_) = A7(^L - /3)^ + As{i,L - P)^ + BritPL - /3) + B8(^l - /3)'. (7.13) 
Moreover, Dl_(f) can be determined by 

DU = ( ^Dl_S + D,_ { I) D,_S] + D,_ ( ^cu) , (7.14) 



7(7-1) V5 ^- \SJ "-J \A{x 

which follows from \7.11\) . The above formulas will be used to resolve the sonic case. See 
Proposition \7.4\ 

We also note that, in order to resolve the contact discontinuity wave Fq, the following 
equation of wq will be used, 

^--i^o^^ifo, B,^( H,^( (7.15) 



dt ■ " dx ^' " \pc' u)' ^ \-^pc'u^ 

Now let us present the LGRPqo solver for problem f l7.6p - (l2.2l) in the following proposition, 
corresponding to the wave configuration in Fig. 17.11 

Proposition 7.2 (Linear GRP). Assume a typical wave configuration for the generalized 
Riemann problem of ^7.6\ ) and h2.2^) as shown in Fig \7.1\ Then {dxQ)*L and {dxQYn ^'^^ 
determined by the set of linear equations 

[(A_/ - B.)VQW_\l{dxQ)l = D,_w4Ul) - H4UI), 

(Ao/ - Bo)UdxWo)l - (Ao/ - Bo)Ud.wo)}, = -Ho{Ul) + Ho{U*^), 

[{VqF - aWQU)iaI - JMO^QTr - D^Ur - Ur) = 

[{VqF - aVQU){aI - J)]l{dxQ)R - [(VqF - aVQt/)//]^^ + [(VqF - aVQU)H]R. 

(7.16) 
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Here, D\_W-{Ul) is determined by \7.9\} , a is the speed of the shock associated with A+ = 
u + c. 

Moreover, for the sonic case where the t-axis is located in the rarefaction, the dtQ at 
t-axis ((3 = 0) are determined by 

[(A_/ - B^)VQW^]{dtQ) = Dx^w^ - H^, 

1 A'{x) (7-17) 

dtu H dtp = . cu, 

pc A[x) 

where the U in the coefficients takes value at t-axis where A_ = 0. 

Proof. The linear system fl7.16p for (9*f/)2 and {d*U)*j^ can be obtained by combining 
(17. 9p . (14. ip and (14. 5 p with / = 1 and using the following expressions 

Dx_ (w.iU)) = {X.I - B.)VQW.d^Q + H_, 

Dao (woiU)) = (Ao/ - Bo)d,Wo + Ho, 

D^[{VqF - aVgU)] = {VqF - aVQU)[iaI - J){d.,Q) + H] - DlaU. 

To resolve the sonic case, we shall use the differential relation along the A_ characteristic 
curve: 

D,_u + 1d._P=^cu, (7.18) 
which is equivalent to (TTTT]) . (ITTTp follows from (EH]) and ( TH^ by setting A_ = 0. □ 

The QGRPoo solver for the nonsonic case and sonic case are presented in the following 
two propositions, respectively. 

Proposition 7.3 (Quadratic GRP: Nonsonic case). Assume a typical wave configura- 
tion for the generalized Riemann problem of jl.O^ and Ii2.2\) as shown in Fig. 7J_. Then 
{'^xQYl ^'^^ i'^xQ)*R ^'^^ determined by the set of linear equations 

Mr{Ul){dlQ)l = Dl_w4Ul) - Br{Ul (8^)1) 

M,(f/2)(9>o)2 - Me(f/^)(9>o)J? = -B,{Ul {8^)1} + B,{U*^, {8^)*^), 

M,{U*^,a){8lQ)l - Dl^m - Ur) = -B,{U*^, {8.,U)l,a,D„a) 

+ M,{Ur, a){8lQ)R + B, {Un, {8,U)r, a, D^a) . 

Here, D\_w.{U^) is determined by a is the speed of the shock associated with A+ = 
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u + c, 

BriU, d.,U) = [Dx_ (A-/ - fi-) - (A_J - B^)d,B^]d^,w^ 

+ (A_/ - B^fd,{VQW^)d,Q + (A_/ - B_)d,H_ + D^_H^, 

M,{U) = (AoJ - 5o)^ 

Bc{U, d^U) = D^iXoI - Bo)d,Wo - (Ao/ - Bo)d,Bod,Wo + (Aq/ - Bo)d,Ho + Z^Ao^^o, 
M,(f/, a) = (VqF - aVQU){aI - J)\ 

B,{U, d,U, a, D^a) = {D^{VqF) - 2D,aVQU - aD,{VQU))D^Q 

+ {VqF - aVQU)[D,{aI - J) - {al - J)d,J]d,Q 
+ {VqF - aVQU)[{aI - J)d,H + D^H]. 

Proof. This proposition can be proved by combining (17.101) . (14.11) and (14. 5 p with / = 2 
and using the following expressions 

Dl {w^{U)) = Mr{U)dlQ + BriU, d,U), 
Dl iwoiU)) = Me(f/)9>o + B,{U, d,U), 

and 

Dl{F{U) - aU) = MsiU, a)dlQ - DjaU + d^U, a, D^a). 

□ 

Proposition 7.4 (Quadratic GRP: Sonic case). Assume that the t-axis is located in- 
side the rarefaction wave associated with A_. Denoting by ^ = {S,(f),ip), then for any point 
= (0, At) with At being small, we have 

$(P.) = $(0) + DxMP*)^t + DL<f (/?*)— + OiAt'), (7.20) 
where (3^ is the root of 

/3 + Da_A_(/3)^ = 0, (7.21) 
and D{ J^, 1=1,2 are determined by ([Op-([7J^. 

The other wave configurations can be treated similarly. In particular, if a A+-rarefaction 
wave is involved, in order to get the linear equation for {d^Q)*^ analogous to the first 
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equation of fl7.19p . it requires one to derive the L(Q)-equation for w^. However, a bet- 
ter choice for us is to use the following property of system (17. 6p : (17. 6p holds true un- 
der the transformation ^ : {p,u,p, A){x,t) — )• {p, —u,p, A){—x,t). In fact, if we de- 
note by g = ^(Q), then Ul = ^(f/*), Ul = ^(?7^), A_(t/^) = A+(f/£). By ex- 
pressing D\_w_{(5) as a function of Ul, {dxU)L, {dlU)L, A'{x), A"{x) and (3: D\_w„{(5) = 
^{Ul, idM)L, idlU)L,A'{x),A"{x),/3), we have 

MriUl -A'{x), A"{x)){dlQ)l =W{Ul, d^Ur., {dlU)^, -A'{x), A"{x), X^UD) 

-Br{Ul{d,U)l-A'{x),A"{x)). (7.22) 

Noting that (d^Q)}^ = ((<9^p)l,, -{d^u)l, (<9^p)2), (I7.22p is indeed the derived linear equation 
for [dlQYj^. The sonic case corresponding to the A+-rarefaction wave can be resolved using 
the same technique. 

7.2. Tests for the GRP solvers 

In this section we assess the performance of the GRP solvers for the compressible Euler 
equations system, i.e. (17. 6p with A{x) = 1. The aim is to show, via several test problems, 
the accuracy and behavior of the present solvers. As tests, we use the generalized Riemann 
problems proposed by [sj and construct new ones with large jumps in pressure. The first test 
has no jump discontinuities in the state variables but admits discontinuities in derivatives 
at the interface. The more demanding test problems are constructed from the first one, 
by adding a discontinuity in pressure. Six new cases are thus generated by varying the 
strength of the initial pressure jump Ap = [p^ — Pr)/pr at the interface, namely Ap = 10'^, 
k = —2, ■ ■ ■ ,3. The last test problem for the sonic case are constructed by adding Am = 28 
to the initial flow velocity of the test case corresponding to Ap = 100. 

In jij , the authors test the first five problems using three type of DRP solvers with only 
partial success. Since no exact solutions are known, the reference solutions are obtained 
numerically, by solving the test problems on very fine mesh on the interval [—1, 1] x [0,to] 
of(x,t). To do this, ^he authors of g , [lH suggest using the Random Choice Method or 
Weighted Average Flux method to avoid the large nonphysical oscillations of early time 
solution. To this aspect, a detailed description can be found in js], [lH, which is beyond the 
scope of this work. Here, our numerical reference solutions are obtained simply by using the 
Godunov flux in the context of finite volume method and then correcting their values on the 
early time interval [0, to/20] using an interpolation method. Such a measure does not affect 
our accuracy tests. 

For each of these tests, we will compare the GRP solver based solution at the interface 
X = 0, as a function of time determined by (12. 5p . with the reference numerical solution. As 
will be shown, the present GRP solvers is truly accurate, having the expected accuracy not 
only for all the test cases in jsf, but also for cases with much larger initial jump in state 
variables. 
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Tabic 7.1: The Ljxj error of U and convergence rate of the acoustic GRP solvers 





t = 0.1 


t = 0.05 


t = 0.025 


t = 0.0125 


Solver 


Error Order 


Error Order 


Error Order 


Error Order 


LGRPi 
QGRPi 


2.420e+0 - 
l.Olle+0 - 


4.407e-l 2.46 
9.861e-2 3.36 


9.592e-2 2.20 
1.127e-2 3.13 


2.251e-02 2.09 
1.439e-3 2.97 




Figure 7.2: Acoustic case: Reference solution and GRP solvers based solutions. 



7.2.1. Acoustic case with continuous state and jump in derivatives 
This test corresponds to the following initial condition 



Pl( 


x,0) 


= 1 + 0.56431X + 2.62896x^ 


ul 


;x,o) 


= 0.03125 - 1.024X + 1.92x2 


Pl{ 


x,0) 


= 10 -0.216X + 1.08x2, 


PR 


:x,o) 


= 1 + 2.04204X, 


Ur 


(x,0) 


= 0.03125 - 0.25X + 0.75x2, 


PR 


>,0) 


= 10, 



(7.23) 



which is indeed the initial conditions used as Test 2 in [8| with slight modification, keeping 
only the same leading terms up to second order at x = 0. The initial condition fl7.23p 
has a continuous state but with discontinuous derivatives at x = 0. The solution for this 
problem contains a left-going and a right-going acoustic waves. The t-axis is located in the 
intermediate region of the two acoustic waves. This test aims at testing the accuracy of 
the acoustic GRP solvers. Fig. 17.21 shows the solution of LGRPi solver and QGRPi solver, 
for each component of U, and the errors measured in L^o with the rate of convergence are 
displayed in Table I7.1[ 

7.2.2. Nonsonic case with jump in initial state 

In this subsection, we test the GRP solvers in the nonsonic case with initial conditions 
having jump in state variables. The initial conditions are generated from fl7.23p by adding 
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Table 7.2: The Loo error of U and convergence rate of the LGRPoo solvers: Nonsonic case 



t = to t = to/2 t = to/4: t = to/8 



Ap 


to 


error Order 


error 


Order 


error 


Order 


error 


Order 


0.01 


0.1 


2.456e+0 - 


4.478e-l 


2.46 


9.762e-2 


2.20 


2.294e-2 


2.09 


0.1 


0.1 


2.782e+0 - 


5.100e-l 


2.45 


1.115e-l 


2.19 


2.630e-2 


2.08 


1 


0.1 


6.012e+0 - 


1.128e+0 


2.41 


2.507e-l 


2.17 


5.881e-2 


2.09 


10 


0.05 


5.823e+0 - 


1.406e+0 


2.05 


3.515e-l 


2.00 


8.778e-2 


2.00 


100 


0.01 


1.265e+l - 


2.810e+0 


2.17 


6.501e-l 


2.11 


1.600e-l 


2.02 


1000 


0.005 


5.201e+2 - 


1.277e+2 


2.03 


3.010e+l 


2.08 


7.300e+0 


2.04 




Table 7.3: 


The Loo error of U and convergence rate of the QGRPoo 


solvers: Nonsonic case 








t = to 


t = to/2 


t = to/4: 


t = to/8 


Ap 


to 


error Order 


error 


Order 


error 


Order 


error 


Order 


0.01 


0.1 


1.024e+0 - 


9.997e-2 


3.32 


1.157e-2 


3.11 


1.517e-3 


2.93 


0.1 


0.1 


1.141e+0 - 


1.113e-l 


3.36 


1.278e-2 


3.12 


1.721e-3 


2.89 


1 


0.1 


2.289e+0 - 


2.250e-l 


3.35 


2.714e-2 


3.05 


3.105e-3 


3.13 


10 


0.05 


1.729e-l - 


2.035e-2 


3.09 


3.411e-3 


2.58 


7.979e-4 


2.10 


100 


0.01 


3.350e+0 - 


4.900e-l 


2.77 


7.000e-2 


2.81 


l.OOOe-2 


2.81 


1000 


0.005 


8.220e+l - 


1.850e+l 


2.15 


2.800e+0 


2.72 


4.000e-l 


2.81 



a term in p^ and thus generating a jump Ap = (Pi(0, 0), Pr(0, 0))/Pr(0, 0) in pressure at 
x = 0. 

The Loo error of vector U with the convergence rate for the LGRPoo and QGRPoo solver 
are tabulated in Tables 17.21 and I7.3[ respectively. We can see that for all cases the LGRPoo 
attains order two and QGPRoo solver is essentially third order. For the QGRP solver, the 
decay of accuracy in some cases may be caused by the limited resolution of reference solution. 
Fig. 17.31 show the results of the acoustic LGRPi (resp. QGRPi) solver in comparison with 
that of the LGRPoo (resp. QGRPoo) solver with Ap arranges from 0.01 to 10. When Ap 
is small, the acoustic solvers do serve as good approximations of their counterpart ones. 
However, as the jump Ap increases, the performance of the acoustic solvers becomes worse. 
As Ap = 10, the acoustic solvers give absolutely wrong initial slopes. Indeed, the behaviors 
of the acoustic solvers are essentially the same with that of the approximate solvers studied 
in [sl, [HI. For even larger pressure jump cases Ap = 100 and Ap = 1000, the solution 
profiles are shown in Fig. 17.41 We can see that, for all these test cases, the LGRPoo and 
QGRPoo solver based solutions agree well with the reference solutions. 

7.2.3. Sonic case 

For the sonic case, the test problem is generated by adding Am = 28 on the initial velocity 
u{x, 0) of the generalized Riemann problem in previous section corresponding to p = 100. 
Compared to the previous tests, it is more difficult to compute the reference solution for this 
case and we need to use a finer mesh with a smaller time interval. The reason is twofold. For 
the first, the solution is singular in the rarefaction fan, and for the second, we have observed 
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Figure 7.3: Nonsonic case: Reference solution and GRP solvers based solutions. Left: LGRPoo and LGRP, 
solver; Right: LGRPi and LGRPi solver. From top to bottom: Ap = 0.01, Ap 0.1, Ap = 1, Ap = 10. 
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Table 7.4: The Loo error of $ and convergence rate of the GRPqo solvers: Sonic case 





t = to 


t = 2/3to 


t = to/2 


t = l/3to 


Solver 


Error Order 


Error Order 


Error Order 


Error Order 


LGRPoo 


1.114e+l - 


4.809e+0 2.06 


2.671e+0 2.06 


1.165 2.05 


QGRPoo 


1.052e+0 - 


3.208e-l 3.13 


1.304e-l 3.13 


4.006e-2 2.91 



an aberration phenomenon when computing reference solution. The aberration phenomenon 
is illustrated by Fig. 17.51 the computed reference of (or E) exhibits a weak discontinuity 
point and an aberration region. However, for the GRI, the 5* and ip, such a phenomenon is not 
observed. This phenomenon is different from the afore- mentioned early-time oscillation jsf, 
since it is GRI-dependent. As the mesh is refined, the weak discontinuous point converges 
to the singularity (0, O"*") and the numerical solution converges. 

This phenomenon can be viewed as a numerical justification of the fact that the second 
time derivative of a variable, expect for the GRI, takes infinite value at the singularity. See 
Section 15.21 

The errors in terms of the vector $ = (5*, ip, 0) and the convergence rates for the GRP 
solvers are displayed in Table 17.41 As suggested in Section 15. 2^ for resolving the sonic case, 
we use the Newton iteration method with initial gauss /3 = to solve f l7.2ip . Here, for the 
tolerance TOL = l.Oe — 7, the number of iterations required for convergence is no more than 
three. 



8. Numerical schemes 

In this section, we turn using GRP solvers to construct one step high order numerical 
schemes, namely, the GRP schemes. In the introduction, we have described the process 
of implementing the second-order numerical scheme, where the LGRP solver provides a 
second-order approximation of the flux function from a piecewise linear discontinuous initial 
data. The process of implementing the QGRP solver based third-order numerical scheme is 
similar. The differences is that we need to provide a third-order subcell data reconstruction 
on each time step and use two point quadrature for the integral of (11. 4p to compute the 
numerical flux, i.e. 

Fj+i/2 = uJiF{U{xj+y2.ri)) + uj2F{xj+y2,r2). (8.24) 

On each quadrature points (xj+i/2,rj), the vector U are calculated through (12. 5p . wherein 
the L'^(0,0"'"), c?ft/(0, 0"*") and d^U{0,0'^) are determined by solving a generalized Riemann 
problem on the cell interface using the QGRP solver. 

In the following, we present several one- dimensional examples to test the performance of 
our schemes. The uniform size meshes are used for all the test cases. For the second-order 



scheme, the van Leer limiter [25| is used to perform the linear reconstruction. For the third- 



order scheme, we use the same reconstruction method as in [18|. In fact, we use the 5rd 
order WENO technique to reconstruct pointwise variables of U at each cell interface. Then 
based on the cell interface values and the cell averages of ?7", a third-order polynomial is 
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Figure 7.5: Sonic case: Reference solution and GRP solvers based solutions. Uniform mesh of 2.0e-4 cell 
size are used for computing reference solution. 
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Exact 
2nd Order 
3rd Order 






- Exact 
2nd Order 
3rd Order 




Figure 8.1: Numerical solutions of Sod problem: 100 grid point are used. 



constructed as the subcell flow distributions at time t". In the following numerical examples, 
the WENO reconstruction is carried out based on the characteristic decomposition 19| and 
the CFL number is set to be 0.5. 

For all the problems, the GRP solutions are plotted against the exact solutions. The 
solid lines represent the exact solution, the circles show the second-order scheme solution, 
while the crosses stand for the third-order scheme solution. 



.1. Sod problem 

The first test is the standard Riemann problem proposed by sod 2l|. The gas is initially 



at rest with p = 1, p = 1 for —5 < x < and p = 0.125, p = 0.1 for < a; < 5. At time 
t = 2, the numerical solutions with 100 points are shown in Fig. 18.11 We can see that both 
of the computed solution agree well with the exact one and the third-order scheme shows 
better performance. 

8.2. 123 problem 

This example was first proposed by j^. The initial data is given with (p, u,p) = (1, 2, 0.4) 
for —5 < X < and {p,u,p) = (1,2,0.4) for < a; < 5. The numerical solutions at time 
t = 1.2 are shown in Fig. 18. 2[ This test case demonstrates the ability of the GRP schemes 
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Figure 8.2: Numerical solutions of 123 problem. 100 grid point are used. 



to preserve the positivity of the density, pressure and internal energy. Again, the internal 
energy profile conforms the better performance of third-order scheme. 

8.3. Woodward- Colella blast wave problem 



This is a problem proposed by [26|. The diatomic gas is initially at rest, and the density 
is unit everywhere. The pressure is p = 1000 for < a; < 10 and p = 100 for 90 < a; < 100, 
while it is only p = 0.01 in 10 < x < 90. Reflecting boundary conditions are applied at 
both ends and the output time is t = 3.8. Numerical solutions with 400 points are shown in 
Fig. l8.3l to exhibit the performance of both schemes. This test case clearly demonstrates the 
capability of both schemes in the capturing of strong shock waves. The third-order scheme 
capture much sharper solution than the second-order scheme in the density and internal 
energy distribution. 

8.4. Shock- density wave interaction 



The Mach 3 shock-entropy wave interaction [20| is specified by the initial condition: 
(p,n,p) = (3.57134,2.629369,10.33333) for < x < 1 and ip,u,p) = (1 + 0.2 sin(A;x), 0, 1) 
for 1 < a; < 10 with k = 5. The solution of this problem consists of a number of shocklets 
and fine scales structures which are located behind a right-going main shock. The computed 
density profile with 400 points, at t = 2.0, is shown in Fig. 18.41 Again the third-order scheme 
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25 50 75 100 25 50 75 100 



Figure 8.3: Numerical solutions of Woodward-Colella blast problem. 400 grid point are used. 
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Figure 8.4: Numerical solutions of shock-density wave problem. 400 grid point are used. 



works better and captures much finer scale structures at high frequency waves behind the 
shock. 

8.5. Steady flow in a converging- diverging nozzle 

We now use the examples in 0, Sect. 6.5] to test the abihty of the GRP schemes to 
attain the steady state of a flow. Consider a flow in a converging- diverging nozzle, which 
occupies the internal < x < 1 and has a smooth cross-sectional area function A{x) given 
by the following expression: 



A{x) 



exp ( - log(Ain) sin^(27rx)) 

Aex exp ( - log(Aex) siu^ ' Mln^l 



< X < 0.25; 
0.25 < X < 1, 



.25) 



where A^-^ = 4.8643 and A^^ = 4.2346. See Fig. 18. 5[ For a steady duct flow of a perfect gas, 
the Mach number M{x) = u{x)/c{x) is determined by A{x) through the algebraic relation 



[A{x)f 



7 + 1 



1 



7 



[M(x)]2 

Then the steady flow profiles in the nozzle are given by 

7-1 



-[M{x)f 



7 + 1 
7-1 



(8.26) 



Pix) = Po 1 + 



p(x) = Po 1 + 



2 

7-1 



[Mix)r 
[M{x)r 



1 

7-1 



7-1 



(8.27) 



u{x) = M{x)^/l'yp{x)/p{x)), 

for the flow being smooth, where po and po need to be specified. 
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Figure 8.5: Nozzle contour. 
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Steady Solution 
2nd Order 
3rd Order 



0.6 - 



0.4 - 



0.2 - 



0.25 



0.5 



0.75 



Figure 8.6: Large time flow in Laval nozzle: Case A. 22 grid points are used, at time t ~ 15.5. 
The initial data we use are 



where pb is a constant determined by the steady state solution at z = 1. We consider two 
cases. In both cases we take Po = Po = ^ and A{x) as in fl8.25l) . 

(A) A smooth flow where p{l) = 0.0272237 is obtained from (18.271) by taking a; = 1 in 
lKm> . leading to M(l) = 3. 

(B) Setting p{l) = 0.4 leads to a discontinuous steady state solution, as shown by sohd 
lines in Fig. 18.81 

We use the strategy in 0, Sect. 6.5] to deal with the boundary conditions at x = and 
1. In both cases, the number of grid points used are 22. As shown in Fig. 18.61 and 18.81 
both of the GRP solutions at t = 15.5 are good agreement with the exact solution. The 
third-order GRP solutions are closer to the analytical solutions than the second-order ones. 
Moreover, as shown in Fig. 18 ■7[ the GRP solutions almost attain the steady state at time 
t = 2.5. This shows that the GRP solutions converges to steady solution quickly. 



U{x,0) 



Ul = (po,0,P6), 
Ur = (poiPb/Po) 



0<x < 0.25, 
V^,0,pfe), 0.25 < a; <1, 
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Figure 8.7: Large time flow in Laval nozzle: Case B. 22 grid points are used, at time t = 2.5. 




Figure 8.8: Large time flow in Laval nozzle: Case B. 22 grid points are used, at time t = 15.5. 
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Appendix A. Formulae in Section 17.11 

For the general case of 7 > 1, the L-equations of w_ = {S, tp) in Section [TTT] yield 

D^_S{0,P) = A,{^L-f3)^, 

Dx_tp{0, /3) = A2{iJL -13)^ + AMitPL - P)^ + ^o(/3)(^L - 
where A2 = A2 — Zq^Pl) and 

+52(^L-/3)^, if 7 7^ 3, 5/3; 
f 4g - /3) + B^ii^L - if 7 = 3; 

B.Mi^L - - 1^ Hi^L - if 7 = 5/3, 



(A.l) 



7+1 



^1 = {i^L-h)~~Dx_S{<d,l3L) 
A2- 



7(37 - 1)^. 



1 (V^L - /3l)^/^a_5(/3l) + (V^L - /3l)"^/^a_^(/3l). 



The coefficients Aj, 5^, are defined as in Table \KA\ and \K72\ 
The function and Z2{P) in Proposition 17. II are as follows 



^2(/3) 



37-1 
7 + 1 



2 



-2A 



17 



+ 



7-1 
27 

Ci 



A 



16 



272^1 



18 



27 



2725L 
2(7-1), 



C4 



7^5^ 

C2 



Ce 

^I^Sl 2 



7+1 2('y — 1) 7-3 



+ 



7 + 1 
7-1 



7-3 



2 

7-1 
2(7 - 2) 



7-3 



2(7-2) 
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Tabic A.l: The coefBcicnts A^, - ■ ■ ,Aig and Bi, 



^3 
A, 

A, 
Ae 
Ar 
As 
A9 
Aio 

An 

A12 

An 
An 

^15 

AlQ 

An 
An 
An 



1 



7(37-1)Sl 



7(7+1)51, 



Ai 
Ai 



^(^3 - A,) 

(3-7)(7+l) A 



27_ ^ 3-7 



, A, + i±^A4) 



7(7-l)5i^l 

1 4 , 

27(7-l)5£ ^7 

As 



27(7-l)5£ 
2 



7+1^ ^ 



2{^^7 



2(7-1)' 
Ag + An 

2A1A12 
2A1A13 

AlAio + A2A13 + ^3^12 
+ v43v4i3 

A2A12 



5l 

^2 
^3 
^4 

B, 
By 

Bs 

B9 

Bn 

Bn 
B12 



7-1 
7-3 



' A{x) 



2(7-1) A'jx) 
(7+l)(37-5) A{x) 
7-1 A'(x) 
7+1 A(x) 

2(7-1) A'(x) 

(7+1)2 

^(51-^3) 



2 — B4) 
_3_75^ - 3±li?3 



■^S2 



2(7-1)- 



2 ^7-1 y 

1 r> 

27(7-1)5l'°7 

1 

27(7-l)Si 



2(7-1)- 



Table A.2: The coefficients Ci, ■ ■ ■ ,6*9 

Ci 2AiB^ 

C2 2AiBn 

C3 A2B, + AnBi-^A-j^ 

C4 A359 + AiBn + Ai3i?i - f Ag^^ 

C, A^Bn + AnB2 -f - f) ^ 

C, A,Bn + Aii?i2 + AnB2 + - ^) 4^ 

^7 BiB,-^^B,^-i.l[^) 

Cs B2B, + B,B,o + (-f 58 + :^ - ^) + Sf^. (i^)' 

^^9 ^2^10 + J-^ ) -A{F) ~ ^ 1^ A(^ j 
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